function answer = s_n(delta, n, k)
answer = vpa(exp(-delta * k), 50);
sum = 0;
for j = 0:n
    sum = vpa(sum + vpa(((-1)^j) * (delta^(n - j)) * (k^(-j - 1)) / factorial(n - j), 50), 50);
end
answer = vpa(answer * sum, 50);
end

